function out = simPerturbation(model,shocks)
% Simulating the following variables (in deviation from steady state)

numSim       = size(shocks,2);
ny           = model.ny;
nx           = model.nx;
nx1          = model.nx1;
nx1          = model.nx1;   %number of endogenous states

ySim         = zeros(ny,numSim);
xSim         = zeros(nx,numSim);

if model.pruningOn == 1
    xf_cu    = zeros(nx,1);
    xs_cu    = zeros(nx,1);
    xrd_cu   = zeros(nx,1);
    xAllSim  = zeros(nx+(model.appOrder-1)*nx1,numSim);
    for s = 1:numSim
        if model.appOrder == 1
            ySim(:,s) = model.g0 + model.gx*xf_cu;
            xSim(:,s) = xf_cu;
            xAllSim(1:nx,s) = xf_cu;
        elseif model.appOrder == 2
            xfxf_cu   = kron3(xf_cu,xf_cu);
            ySim(:,s) = model.g0 + model.gx*(xf_cu+xs_cu) + model.Gxx*xfxf_cu + 0.5*model.gss;
            xSim(:,s) = xf_cu+xs_cu;
            xAllSim(1:nx+nx1,s) = [xf_cu;xs_cu(1:nx1,1)];
        elseif model.appOrder == 3
            xfxf_cu   = kron3(xf_cu,xf_cu);
            xfxs_cu   = kron3(xf_cu,xs_cu);
            ySim(:,s) = model.g0 + model.gx*(xf_cu+xs_cu+xrd_cu) + model.Gxx*(xfxf_cu + 2*xfxs_cu) ...
                        + model.Gxxx*kron3(xf_cu,xfxf_cu) ...
                        + 0.5*model.gss + 3/6*model.gssx*xf_cu;
            xSim(:,s) = xf_cu+xs_cu+xrd_cu;
            xAllSim(1:nx+2*nx1,s) = [xf_cu;xs_cu(1:nx1,1);xrd_cu(1:nx1,1)];
        end
        
        % Updating the states (dev from steay state)
        xf_cup = model.hx*xf_cu + model.eta*shocks(:,s);
                         
        xf_cu = xf_cup;
        if model.appOrder > 1
            xs_cup = model.hx*xs_cu + model.Hxx*xfxf_cu + 0.5*model.hss;
            xs_cu  = xs_cup;
        end
        if model.appOrder > 2
            xrd_cup = model.hx*xrd_cu + model.Hxx*(2*xfxs_cu) ...
                + model.Hxxx*kron3(xf_cu,xfxf_cu) ... 
                + 0.5*model.hss + 3/6*model.hssx*xf_cu;
            xrd_cu  = xrd_cup;
        end 
    end
else
    x_cu = zeros(nx,1);
    for s = 1:numSim
        xSim(:,s)   = x_cu;
        if model.appOrder == 1
            ySim(:,s) = model.g0 + model.gx*x_cu;
            x_cup     = model.hx*x_cu + model.eta*shocks(:,s);
        elseif model.appOrder == 2
            xx_cu     = kron3(x_cu,x_cu);
            ySim(:,s) = model.g0 + model.gx*x_cu + model.Gxx*xx_cu + 0.5*model.gss;
            x_cup     = model.hx*x_cu + model.Hxx*xx_cu + 0.5*model.hss + model.eta*shocks(:,s);
        elseif model.appOrder == 3
            xx_cu     = kron3(x_cu,x_cu);
            ySim(:,s) = model.g0 + model.gx*x_cu + model.Gxx*xx_cu + model.Gxxx*kron3(xx_cu,x_cu) ...
                        + 0.5*model.gss + 3/6*model.gssx*x_cu;
            x_cup     = model.hx*x_cu + model.Hxx*xx_cu + model.Hxxx*kron3(x_cu,xx_cu)...
                        + 0.5*model.hss + 3/6*model.hssx*x_cu + model.eta*shocks(:,s);
        end
        % Updating the states       
        x_cu  = x_cup;
    end
    
end

%% Saving output
out.y_cu  = ySim;
out.x_cu  = xSim;
if model.pruningOn == 1
    out.xAll_cu = xAllSim;
end

end